Charged Ising Model of Neutron Star Matter 
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Background: The inner crust of a neutron star is believed to consist of Coulomb-frustrated complex structures 
known as "nuclear pasta" that display interesting and unique low-energy dynamics. 

Purpose: To elucidate the structure and composition of the neutron-star crust as a function of temperature, 
density, and proton fraction. 

Methods: A new lattice-gas model, the "Charged-Ising Model" (CIM), is introduced to simulate the behavior 
of neutron-star matter. Preliminary Monte Carlo simulations on 30 3 lattices are performed for a variety of 
temperatures, densities, and proton fractions. 

Results: Results are obtained for the heat capacity, pair-correlation function, and static structure factor for a 
variety of conditions appropriate to the inner stellar crust. 

Conclusions: Although relatively simple, the CIM captures the essence of Coulomb frustration that is required 
to simulate the subtle dynamics of the inner stellar crust. Moreover, the computationally demanding long-range 
Coulomb interactions have been pre-computed at the appropriate lattice sites prior to the start of the simulation 
resulting in enormous computational gains. This work demonstrates the feasibility of future CIM simulations 
involving a large number of particles as a function of density, temperature, and proton fraction. 

PACS numbers: 05.50. +q, 26.60.-c, 26.60.Gj, 26.60.Kp, 51.30.+i, 95.30.Tg 
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I. INTRODUCTION 

Neutron stars are compact objects with radii of the 
order of ten kilometers and masses comparable to that 
of the Sun. The solution of the Tolman-Oppcnheimer- 
Volkoff (TOV) equations [l|, j which prescribes the 
structure of spherically symmetric, self-gravitating com- 
pacts object in hydrostatic equilibrium, provides infor- 
mation on the density profile of the star. Remarkably, 
the structure of neutron stars depends exclusively on the 
nuclear equation of state (EOS). Given the constraint 
of hydrostatic equilibrium, the density profile spans an 
enormous range of densities: from the extremely dilute 
crustal densities up to core densities that may greatly 
exceed nuclear-matter saturation density. Understand- 
ing what novels phases of matter emerge under these ex- 
treme conditions is both fascinating and unknown 0, . 
Moreover, it represents one of the grand challenges in 
nuclear physics: "How docs subatomic matter organize 
itself?"^. 

The highest density attained in the stellar core depends 
critically on the equation of state of neutron-rich matter. 
Although at such high densities the EOS is poorly con- 
strained, it has been speculated that many exotic phases 
may emerge under such extreme conditions. These may 
include pion or kaon condensates 0, 01, strange quark 
matter Q , and color superconductors [HJ ■ It is also 
often assumed that the uniform core may have a non- 
exotic component consisting of neutrons, protons, elec- 
trons, and muons in chemical equilibrium. However, at 
densities of about half of nuclear-matter saturation den- 
sity, the uniform core becomes unstable against cluster 
formation. At these "low" densities the average inter- 
nucleon separation increases to such an extent that it 
becomes energetically favorable for the system to segre- 



gate into regions of normal density (nuclear clusters) and 
regions of low density (neutron vapor). The transition 
region between the homogeneous and non-homogeneous 
phases constitutes the crust-core interface. It is the aim 
of this work to study the structure and composition of the 
crust-core interface where distance scales are such that 
the Coulomb and nuclear interactions become compara- 
ble in strength. Under these unique conditions neutron- 
rich matter becomes "frustrated" . Frustration, a preva- 
lent phenomenon characterized by the existence of a very 
large number of low-energy configurations, emerges from 
the impossibility to simultaneously minimize all elemen- 
tary interactions in the system. In the inner stellar crust 
this leads to a myriad of complex structures — collectively 
known as "nuclear pasta" — that are radically different in 
topology yet extremely close in energy. Moreover, due to 
the preponderance of low-energy states, frustrated sys- 
tems display an interesting and unique low-energy dy- 
namics. For example, it has been speculated that pasta 
formation could enhance the coherent scattering of neu- 
trinos from such exotic structures. This could have im- 
portant consequences on the supernova explosion mech- 
anism and subsequent cooling dynamics [lll - ll3| . 

In this contribution we are interested in the equation 
of state of neutron-rich matter at densities of relevance 
to the inner stellar crust [3, EH- We will model this 
charge-neutral system in terms of its basic constituents, 
namely, neutrons, protons, and an ultra-relativistic Fermi 
gas of electrons. In particular, no ad-hoc biases will be 
introduced in regard to the structure of the exotic pasta 
shapes (i.e., whether they form droplets, rods, slabs, bub- 
bles, etc.). Rather, we will allow the clustering to develop 
dynamically from an initial (random) configurations of 
nucleons. The aim of this work is to explore the dynamics 
of the system as a function of density, temperature, and 
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proton fraction. Note that the original work by Ravenhall 
and collaborators was carried out at zero temperature in 
a mean-field approach fl6l - [l8j . Recently, more sophisti- 
cated approaches — based on Monte Carlo and Molecu- 
lar Dynamics simulations [Tl[ [HI, [l9l - l29| , a Dynamical- 
Wavelet approach [30l - [3^ | . relativistic mean- field calcu- 
lations [33l - l37l | . and Skyrme-Hartree-Fock methods ;38| — 
have been implemented and have confirmed the existence 
of these exotic phases at very low temperatures and mod- 
erate proton fractions. However, given that chemical 
equilibrium suggests that the proton fraction in the inner 
stellar crust is very low — indeed, significantly lower than 
normally assumed — it has recently been put into ques- 
tion whether pasta formation is even possible in such 
proton-poor environments (39l| . Moreover, simulations at 
different temperatures are both critical and interesting 
because the long-range Coulomb interaction is respon- 
sible for the extreme fragility of crystals. That is, the 
melting (or charge-ordering) temperature in crystals T c 
is significantly smaller than the relevant Coulomb energy 
scale -Ecoui = e 2 /a (here a is the lattice spacing). Such 
an energy mismatch introduces a large temperature gap 
(kBT c <kBT^Ecoui) where the system displays uncon- 
ventional pasta-like behavior that reflects the strong frus- 
tration induced by the long-range interactions. In par- 
ticular, condensed-matter simulations with long-range in- 
teractions have reported the opening of a pseudogap in 
the density of states in response to the strong frustra- 
tion (4(j| . This unconventional pseudogap region medi- 
ates the transition from the Wigner Crystal to the Fermi 
liquid. Interestingly enough, the pseudogap disappears 
for a system with only short-range interactions. 

As an alternative to the numerically intensive Monte 
Carlo and Molecular Dynamics simulations, we introduce 
here the "Charged Ising Model" (CIM). The CIM is a 
lattice-gas model that while simple in its assumptions, 
retains the essence of Coulomb frustration. Numerical 
simulations based on this model are not as computation- 
ally demanding because the long-range Coulomb inter- 
action, computed here exactly via an Ewald summation, 
may be pre-computed at the appropriate lattice sites and 
then stored in memory prior to the start of the simula- 
tion. This represents an enormous advantage when try- 
ing to simulate systems with a large number of particles 
as a function of temperature, density, and proton frac- 
tion. It is an important goal of this work to extend earlier 
(fixed-temperature) approaches by studying the thermal 
properties of the crust. Specifically, we rely on classical 
Monte Carlo simulations of the CIM to investigate phase 
transitions in stellar matter in the presence of Coulomb 
frustration. The CIM is reminiscent of an earlier ap- 
proach developed in Refs. (4lTj43j . Yet, it improves on it 
in two respects: (a) by including explicitly the isospin 
degree of freedom that is required for a proper treatment 
of asymmetric matter and (b) by using the Ewald sum- 
mation to properly treat the long-range Coulomb inter- 
action. Although limited in their treatment of quantum 
fluctuations, classical simulations like the ones proposed 



here are essential to uncover correlations that go beyond 
mean-field approaches. In particular, both spatial and 
thermal correlations — as embodied in the static structure 
factor and heat capacity — will be computed as a function 
of density, temperature, and proton fraction in the search 
of signatures of phase transitions. 

The paper has been organized as follow. In Sec.HH the 
CIM and the general framework will be introduced. Re- 
sults from the simulations will be presented in Scc. lIIII 
and then compared against previous findings reported in 
Ref. [3^]. Moreover, we will extend this earlier work by 
following the evolution of the pasta structures as a func- 
tion of the temperature at fixed density and proton frac- 
tion. Finally, we offer our conclusions and suggestions 
for future work in Sec. lIVI 

II. THE CIM MODEL: GENERAL 
FRAMEWORK 

The main constituents of the stellar crust arc neutrons, 
protons, and a background gas of neutralizing electrons. 
At the densities of relevance to the inner crust, the elec- 
trons may be treated as as an ultra-rclativistic Fermi gas, 
namely, with a dispersion relation e(p) = p. Although 
the CIM presented here represents a simplification of 
the model first introduced in Ref. [ll|, it still retains the 
essence of Coulomb frustration, namely, competing in- 
teractions consisting of a short-range nuclear interaction 
and a long-range Coulomb potential. The CIM assumes 
that nucleons are allowed to occupy only the discrete sites 
of a three-dimensional cubic lattice of volume V = L 3 con- 
taining a total number of S sites. The electrons on the 
other hand arc assumed to provide a uniform neutralizing 
background. 

The potential energy consists of a sum of a short-range 
interaction between nucleons and a long-range Coulomb 
interaction between protons and the uniform electron 
background. That is, 

^Total = ^Nuclear + ^Coulomb • (1) 

For the short-range nuclear interaction the potential en- 
ergy is assumed to be given by a sum of two-body terms 
that act exclusively over nearest neighbors. That is, 

1 S 

^Nuclear = 9 Vi 3 nin 3 ' ( 2 ) 

<i,J> 

where m = 0, 1 denotes the occupation number of site i 
and the "elementary" two-body interaction is given by 

Vij = (b + CTiTj) . (3) 

Here Ti is the isospin of the nucleon occupying site i, with 
Ti = +1 for protons and n=—l for neutrons. Note that 
the repulsive short-range nature of the NN interaction 
is simulated here by precluding the double occupancy of 



3 



lattice sites. Also note that the two-body interaction is 
assumed to be isospin dependent to simulate quantum 
statistics. For example, in order to prevent pure neutron 
matter to be bound, the neutron-neutron interaction has 
to be made repulsive, namely, v nn = (b + c) > 0. Indeed, 
we now describe the procedure employed to fixed the two 
parameters b and c. We assume that a completely filled 
lattice containing A — S nuclcons corresponds to nuclear- 
matter at saturation density. That is, 
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a = 1.842 fm. 



(4) 



For such a filled lattice the energy per nucleon of symmet- 
ric nuclear matter and pure neutron matter at saturation 
density are given by 



EsNM 

A 

EpNM 



= 3(6- c) = -16.5 MeV, (5a) 
= 3(6 + c) = +15.5 McV. (5b) 



This choice fixes the two model parameters to the fol- 
lowing values: 



— MeV = -0.167 MeV. 
6 

-^MeV = 5.333 MeV, 

- > 



or equivalcntly, 



v pn = v np = —5.500 MeV. 



= v nn = +5.167 McV. 



(6a) 
(6b) 



(7a) 
(7b) 



To illustrate the dynamics behind this very simple 
choice wc display in Fig JT] results for the energy per nu- 
cleon of infinite nuclear matter (with the Coulomb in- 
teraction turned off and no electrons) as a function of 
both the filling fraction p/p =A/S and the proton frac- 
tion x p = Z/A. Monte-Carlo simulations were performed 
on a cubic lattice of S = (20) 3 sites and at a tempera- 
ture of T w 0. Note that all simulations were started at 
the high temperature of T — 20 McV and slowly cooled 
down to T ~ until the configuration was frozen. The 
EOS for symmetric nuclear matter (x p = 0.5) yields, by 
construction, a binding energy per nucleon at saturation 
density of 16.5 MeV and decreases to about 8 MeV at 
very low densities — corresponding to the binding energy 
of an isolated (symmetric) cluster. Note that in contrast 
to mean-field descriptions that assume nuclear matter 
to be uniform — and thus the energy to vanish at very 
low densities — the lattice-gas model takes full account (at 
least classically) of clustering correlations. At the other 
extreme (x p = 0) pure neutron matter is unbound at all 
densities and yields, by construction, an energy per neu- 
tron at saturation density of 15.5 MeV; this corresponds 
to a symmetry energy at saturation density of 32 McV. 
Note that in the lattice model the energy of pure neutron 
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FIG. 1. (color online) Energy per nucleon as a function of 
filling fraction p/p —A/S and proton fraction x p — Z/A for 
infinite nuclear matter. A proton fraction of x p — 0.5 repre- 
sents symmetric nuclear matter whereas x p — corresponds 
to pure neutron matter. Monte Carlo simulations were per- 
formed on a lattice with S = (20) 3 sites as the temperature 
T^0. 



matter vanishes at half filling (and below) as the lowest 
energy configuration consists of neutrons surrounded by 
empty sites. 

For the emergence of frustration and the concomi- 
tant development of pasta structures, the Coulomb re- 
pulsion between protons is of critical importance. As 
mentioned earlier, at densities of relevance to the bot- 
tom layers of the inner crust (i.e., 10 13 -10 14 g/cm 3 ) the 
competition between the short-range nuclear attraction 
and the long-range Coulomb repulsion is the main driving 
force behind frustration. Whereas in earlier publications 
we adop ted an approximate screened Coulomb interac- 
tion [llf, more recently [39[ we have treated the problem 
exactly by means of an Ewald summation [44| ■ We follow 
the exact Ewald treatment here as well. 

Using Ewald's method we can cast the Coulomb po- 
tential as a sum of two-body interactions plus a constant 
term. That is, 



V c 



oulomb 



= Vn 



p p 



(8) 



where 



,p — 



nj(l+Tj)/2 denotes the proton occupation 



number of site i, v = e 2 /L sets the Coulomb energy scale, 
and Vo is an overall constant j39j j . The dimensionless two- 
body potential uij may be written in terms of short- and 
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long-range contributions: 



form: 



evfc(sij/s ) 



¥0 



exp(— 7r 2 So I 2 ) 



CXp( — 2-Kl\ ■ Sij) , 



(9) 



where 1 = {l x ,ly^z) represents a triplet of integers and 
is the separation between lattice sites i and j in 
dimensionless units. We now proceed with a brief ex- 
planation of the various terms; for a more detailed ac- 
count see Ref. (39|. The Coulomb potential is an inter- 
action with no intrinsic scale. Ewald introduced a scale 
into the problem by adding Z positive and Z negative 
smeared charges at the exact location of each proton. It 
is both customary and convenient to introduce a gaus- 
sian charge distribution with a smearing parameter a s ; 
in the above expression sq = ci s /L. The role of each neg- 
ative charge is to fully screen the corresponding point 
proton charge over distances of the order of the smearing 
parameter. Thus, as long as a s is significantly smaller 
than the box length L, the resulting (screened) two-body 
potential [erfc(s/so)/s] will become short ranged and thus 
amenable to be treated using the minimum-image con- 
vention (45), |4(| . What remains then is a periodic sys- 
tem of smeared positive charges together with the neu- 
tralizing electron background. Whereas in configuration 
space this long-range contribution is slowly convergent, 
the great merit of the Ewald construction is that it can 
be made to converge rapidly if evaluated in momentum 
space, namely, as a Fourier sum. Indeed, the Fourier sum 
is rapidly convergent because (dimensionless) momenta 1 
satisfying Isq ^ 1 make a negligible contribution to the 
Fourier sum. Hence, by suitably tuning the value of the 
smearing parameter, the evaluation of the Coulomb po- 
tential may be written in terms of two rapidly convergent 
sums; one in configuration space and one in momentum 
space [1^] . This is the enormous advantage of the Ewald 
construction. Another enormous advantage — now spe- 
cific to the lattice model — is that one may pre-compute 
the two-body Coulomb interaction Uij for all different 
pairs of lattice sites and then stored them in an array for 
later retrieval during the simulation. 

In what follows we employ a canonical ensemble to 
perform numerical simulations of a system consisting of 
A nuclcons, Z = x p A, protons, and temperature T. A 
configuration in the system may be specified by a col- 
lection of S occupation numbers a = (ai, a%, . . . , as), 
where at each site ai = {p, n, 0}, depending on whether 
the site is occupied by either a proton or a neutron, or 
it remains vacant. Given that the potential energy is in- 
dependent of momentum, the partition function for the 
system factors into a product of a partition function in 
momentum space — that has no impact in the computa- 
tion of momentum-independent obscrvablcs — times a co- 
ordinate space (or interaction) partition function of the 



Z(A, x p , T) = cxp (-/3 VWi(a) 



(10) 



where /3= (fc^T) -1 is the inverse temperature. In turn, 
the expectation value of any momentum-dependent ob- 
servable O may be estimated by performing the appro- 
priate statistical average. That is, 



{O) =Y J 0{a)P a (T) , 



(11) 



where P a (T) represents the probability of finding the sys- 
tem in a given configuration a. In the canonical ensemble 
such a probability is proportional to the properly normal- 
ized Boltzmann factor: 



exp 



-/?Vrotal(<*)) 



Z(A,x p ,T) 



(12) 



Given that the momentum-independent interactions have 
no impact on the kinetic energy of the system, the ex- 
pectation value of the kinetic energy reduces to a sum of 
a classical contribution for the nucleons and a quantum 
contribution for the electrons. That is, 



(K) = ^-Ak B T + -Zk F 
2 4 



2t£ 
3 



T 
% 



(13) 



where k^ = fee^F is the electronic Fermi momentum. The 
total energy of the system is then given by 

(E(A,x p ,T)) = (K(A,x p ,T)) + {Vi om {A,x p ,T)) . (14) 

We note that the sum over a in Eq. ([TTj) runs over a total 
number of configurations given by 



C(A,Z) 



S\ 



Z\{A-Z)\{S-A)\ 



(15) 



This number becomes astronomical even for systems of 
moderate size. Thus, to properly sample the statistical 
ensemble, we rely on a Metropolis Monte-Carlo algorithm 
H3 to generate configurations distributed according to 
Eq. (fT2|) . Given that the kinetic energy of the system 
corresponds to that a classical ideal gas of nucleons and 
an ultra-relativistic Fermi gas of electrons, their contri- 
bution to the heat capacity is both known and smooth. 
Thus, any non-analytic behavior associated with the ex- 
istence of a phase transition must arise from the interac- 
tions. For example, in the case of the heat capacity the 
potential energy contribution will be estimated from the 
fluctuations in the potential energy. That is, 



n 2 Z 



Total) 



(V- 



Total/ 



(16) 



Whereas the heat capacity accounts for the mean- 
square energy fluctuations — which diverge near phase 
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transitions — the static structure factor S(k) provides 
a complimentary observable associated with the mean- 
square density fluctuations [H)]. Moreover, 5(k) is in- 
timately related to a quantity particularly suitable to 
be modeled in computer simulations, namely, the pair- 
correlation function g(r). Indeed, S(k) and g(r) are 
simply Fourier transforms of each other. The pair- 
correlation function g(r) is particularly simple to sim- 
ulate as it represents the probability of finding a pair of 
particles separated by a fixed distance r. For a system 
containing N particles and confined to a simulation vol- 
ume V, g(r) may be computed exclusively in terms of the 
instantaneous positions of the particles. That is, 



7(r) = 1 



V 



N(N - 1) 



i¥=3 



(17) 



where =r^ — r 3 and the "brackets" represent an en- 
semble average. Note that g(r) is normalized to 1 at very 
large distances. Whereas for a uniform fluid the one-body 
density is constant, interesting two-body correlations 
emerge as a consequence of interactions. For example, 
the characteristic short-range repulsion of the TV TV inter- 
action precludes particles from approaching each other. 
This results in a pair-correlation function that vanishes 
at short separations. In the particular case of the CIM, 
this short-range repulsion is enforced by precluding two 
nucleons from occupying the same site. The static struc- 
ture factor is obtained from the pair-correlation function 
through a Fourier transform. That is (49j . 



S(k) = l + y J d\(g{v)-l 



(18) 



Given that the static structure factor accounts for the 
mean-square density fluctuations in the ground state, 
it becomes a particularly useful indicator of the crit- 
ical behavior associated with phase transitions — which 
themselves arc characterized by the development of 
large (i.e., macroscopic) fluctuations. Indeed, the spec- 
tacular phenomenon of "critical opalescence" in fluids 
is the macroscopic manifestation of abnormally large 
density fluctuations — and thus abnormally large light 
scattering — near a phase transition [50j . In this regard, 
the static structure factor at zero-momentum transfer 
provides a unique connection to the thermodynamics of 
the system [5(|. That is, 



S(k = 0) = 



(N 2 ) - (N) 2 _ (N)k B T 



(N) 



V 



(19) 



where kt is the isothermal compressibility of the system. 
The isothermal compressibility is reminiscent of the heat 
capacity [Eq. (|16[)] that accounts for energy rather than 
density fluctuations. As such, they both play a critical 
role in identifying the onset of phase transitions. 

As mentioned earlier, the various configurations of the 
system will be generated via a Metropolis Monte-Carlo 
algorithm with a weighting factor determined by the total 



potential energy of the CIM [see Eq . (TL21 1 . The Metropo- 
lis algorithm is very well known [471 . |49J , so we only pro- 
vide a brief review of those parts of relevance to our im- 
plementation. In particular, all Monte-Carlo moves must 
be consistent with the specified baryon number A and 
proton fraction x p . Thus, given a current configuration 
a, we propose a move to a new configuration a' by se- 
lecting two lattice sites (i and j) at random and then 
simply exchange their occupancies (i.e., a^-f^oy). This 
move ensures that both the baryon number the proton 
fraction are conserved during the simulation. The new 
configuration is accepted provided 

^gUrand, (20) 

where "rand" is a random number between and 1 drawn 
from a uniform distribution. Otherwise, the move is re- 
jected and the original configuration a is kept. 

Initially, the lattice is populated by placing Z = x p A 
protons and N = A — Z neutrons at random throughout 
the S lattice sites. Given that each lattice site is occu- 
pied by at most one nucleon, a total of S— A sites remain 
empty. The simulation starts by thermalizing the sys- 
tem at a temperature that is significantly higher than 
the target temperature T; this prescription prevents the 
system from getting trapped in a local minimum. Once 
the system is properly thermalized at the higher temper- 
ature, a very slow cooling schedule is enforced until the 
desired temperature T is reached. Note that without a 
proper cooling schedule, a system that should crystallize 
at low temperature may end up resembling an amorphous 
solid. Once the system reaches the target temperature T, 
one proceeds to accumulate statistics in order to compute 
the thermal averages for a variety of physical observables. 
However, a strong correlation is likely to exist between 
two neighboring configurations a and a' since they differ 
by (at most) the permutation of two occupation num- 
bers. This correlations can significantly bias the results 
and may lead to an improper estimate of the Monte Carlo 
errors. To prevent this situation from developing, one 
selects uncorrected events by calculating the normalized 
auto-correlation function of a suitable observable O. For 
a large sequence of configurations {a lf a2, . . .}, the auto- 
correlation function of O is defined by the following ex- 
pression: 



fo( m ) 



£„ =1 (On-(0)) (o n+m -(o) 
E„=i (0 n -(0)) 2 



(21) 



where O n = 0(a n ). The decorrelation "time" t is de- 
fined by the condition / (r) =0.1. In Fig. [2] we display 
the auto-correlation function for the total potential en- 
ergy Vrotai for a filling fraction A/S = 0.2, a proton frac- 



tion of 



0.3, and temperatures of T = lOMeV and 



T=15MeV. At the lower temperature it becomes more 
difficult to explore the full energy landscape, thereby re- 
sulting a in a longer decorrelation time. For this partic- 
ular case, T, t . = 7, 218 and t,. = 5, 550. In what follows, 
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all our results are reported with a proper treatment of 
Monte Carlo errors. 
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FIG. 2. The auto-correlation function for the total potential 
energy at fixed density, fixed proton fraction, and two values 
of the temperature. The decorrelation time r is defined by 
the condition /(r) = 0.1. 



III. RESULTS 

We start this section by providing a baseline CIM cal- 
culation that aims to reproduce the results reported re- 
cently in Ref. [3{|. Recall that in Ref. [39[ the temper- 
ature was fixed at T = 1 MeV in order to simulate the 
quantum zero-point motion. It is the goal of our present 
lattice-gas simulation to improve on such a work by ex- 
amining the role of the temperature on the structure and 
dynamics of the inner stellar crust. Ultimately then, this 
sort of simulations will help us explore the phase dia- 
gram as a function of temperature, density, and proton 
fraction. 

Given that the static structure factor at zero mo- 
mentum transfer accounts for density fluctuations [see 
Eq. (|19[) ]. we begin this section by displaying in Figs. [3} 
|5] pair correlations functions and static structure fac- 
tors for neutrons and protons at a fixed temperature of 
T = 1 MeV. Note that because of the discrete nature of 
the lattice, all distances between sites are "quantized" . 
Moreover, due to the periodicity of the lattice, the al- 
lowed values of the momenta are given as follows: 

k=^l (k =0,1,..., 5,-1) , (22) 



Results are presented as a function of the proton frac- 
tion for a lattice of S= (30) 3 sites, and a filling fraction 
of p/p Q =A/S= 0.1875 (or A~ 5,000 nucleons). The pair 
correlation function is characterized by a set of discrete 
peaks at the allowed distances on the lattice. For exam- 
ple, at this relatively low filling fraction, the dynamics 
favors the formation of neutron-rich clusters immersed 
in a dilute neutron vapor (see Fig. lTTl) . Given that the 
neutron-proton interaction is attractive, nucleons orga- 
nize themselves within a cluster by occupying alternating 
lattice sites. Thus, the closest distance between nucle- 
ons of the same species is r m ' m = \p2~a = 2.605 fm, where 
a = 1.842 fm is the lattice spacing [see Eq. ©]. The 
largest peak in both Figs.[3Jand|l]reflcct this behavior. In 
the case of protons — where no dilute vapor is formed — 
other peaks corresponding to more distant protons are 
clearly discernible at distances of 2a = 3.684, \/6a = 4.512, 
yf%a = 5.210, \/l0a = 5.825, ... In the case of neutrons, 
the existence of a dilute neutron vapor gives rise to addi- 
tional peaks and to significant pair-correlation strength 
at larger distances. The corresponding static struc- 
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FIG. 3. Proton pair correlation function for different proton 
fractions at a fixed density of p/p = 0.1875 and a temperature 
of T=lMeV. 

ture factors for both protons and neutrons are shown in 
Figs. [5] and [H respectively. Our results reproduce quali- 
tatively those reported in Ref. [39|. That is, S(k) displays 
a prominent peak that becomes progressively higher with 
increasing proton fraction. The peak in S(k) occurs at a 
momentum transfer k for which the probe (e.g., electrons 
in the case of protons and neutrinos in the case of neu- 
trons) can most efficiently scatter from the density fluc- 
tuations in the system. In particular, if the wavelength of 
the probe is large as compared with the size of the pasta 
structures, the scattering may be coherent. This can 



7 



i 



X D =0.25 

xj=0.2 

X^=0.15 

Xj=0.1 

Xp=0.05 



p/p =0.1875 
T=1MeV 



lid ^fcjj^yLutoiiii 



— 



8 10 12 14 16 18 20 
r[fm] 



CO 




X D =0.25 

Xj=0.2 

Xq=0.15 

xj=0.1 

Xp=0.05 



p/p =0.1875 
T=1MeV 



8 10 12 14 16 18 20 
l=kl_/27C 



FIG. 4. Neutron pair correlation function for different proton 
fractions at a fixed density of p/p =0.1875 and a temperature 
of T=lMeV. 



FIG. 6. Neutron static structure factor for different proton 
fractions at a fixed density of p/p a = 0.1875 and a temperature 
of T=lMeV. 



significantly enhance the response or, equivalently, sig- 
nificantly reduce the electron/neutrino mean-free path. 
Finally, as in Ref . (39j , there is no visible enhancement in 
S(k) at zero momentum transfer as would be expected 
from the putative phase transition from a Wigner Crys- 
tal to a pasta phase. Although it is gratifying that 
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FIG. 5. Proton static structure factor for different proton 
fractions at a fixed density of p/p — 0.1875 and a temperature 
of T=lMeV. 



the CIM reproduces the trends reported in Ref. [39j . an 
important goal of the present work is to explore the evo- 
lution of the system — particularly the dissolution of the 
pasta — as a function of temperature. In analogy to the 
static structure factor that captures the density fluctua- 
tions in the system, we now examine thermal fluctuations 
through a study of the heat capacity [see Eq. (fT6|) ]. We 
start by displaying in Figs.lTTITUl the neutron and proton 
pair correlation functions and corresponding static struc- 
ture factors at a fixed density of pj p a =0.2, a fixed proton 
fraction of x p = 0.3, and for temperatures ranging from 
T=lMeV to T = 6MeV. In particular, note that under 
these conditions of density and proton fraction — and at 
the low temperature of T = 1 MeV — the existence of a 
pasta phase has been well established [ll|, H3 i see a l so 
Fig.EJ 

The T = 1 MeV results set the baseline as these can be 
directly compared against our earlier findings. As before, 
at low temperatures (T < 3 MeV) the system displays the 
strong clustering correlations characteristic of the pasta 
phase. However, as the temperature increases and the 
thermal energy becomes comparable to the binding en- 
ergy per nucleon of the neutron-rich clusters the behavior 
changes dramatically. The large peaks in both the pair 
correlation function and the static structure factor get 
significantly reduced as the system reaches a tempera- 
ture of T~4MeV and both become essentially structur- 
less at T > 5 MeV. Note also the appearance of a small 
peak in g n (r) at r = a as the entropic contribution starts 
to become as important, if not more, than the energy 
contribution. To illustrate the behavior of the system 
as a function of temperature we display in Figs. llUfTBI 
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FIG. 7. Proton pair correlation function for various temper- 
atures, at a fixed density of p/p =0.2, and a proton fraction 



of x p = 0.3. 



FIG. 9. Proton static structure form factor for various tem- 
peratures, at a fixed density of p/p =0.2, and a proton frac- 
tion of a;„ = 0.3. 
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FIG. 8. Neutron pair correlation function for various temper- 
atures, at a fixed density of p/p =0.2, and a proton fraction 



of a; p = 0.3. 



FIG. 10. Neutron static structure form factor for various tem- 
peratures, at a fixed density of p/p =0.2, and a proton frac- 
tion of x p = 0.3. 



Monte-Carlo snapshots at a density of p/p =0.2, a pro- 
ton fraction of x v = 0.3, and temperatures of T— 1 MeV, 
T = 3 MeV, and T = 18 MeV. One can see the gradual 
transition in the structure of the system. At T = 1 MeV 
the system displays the existence of neutron-rich clusters 
surrounded by a dilute neutron vapor. As the tempera- 



ture is increased to T= 3 MeV, some of the weakly bound 
neutrons in the clusters join the vapor and one sees a co- 
existence between the clusters and the vapor. Finally, at 
the very large temperature of T = 18MeV no spatial cor- 
relations remain as the system has been fully vaporized 
into a classical gas of nucleons. Quantitatively, this be- 
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FIG. 11. Monte-Carlo snapshot of neutron-star matter at a 
temperature of T = 1 MeV, a density of p/p = 0.2, and a 
proton fraction of x p = 0.3. The blue and red dots are used 
to display the location of neutrons and protons, respectively. 



FIG. 13. Monte-Carlo snapshot of neutron-star matter at a 
temperature of T = 18 MeV, a density of p/p = 0.2, and a 
proton fraction of x p = 0.3. The blue and red dots are used 
to display the location of neutrons and protons, respectively. 
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FIG. 12. Monte-Carlo snapshot of neutron-star matter at a 
temperature of T = 3 MeV, a density of p/p Q = 0.2, and a 
proton fraction of x p = 0.3. The blue and red dots are used 
to display the location of neutrons and protons, respectively. 



havior is captured by the heat capacity which has been 
computed by calculating the fluctuations in energy [see 
Eq. (|16p ] and is displayed in Fig.QjJ The energy fluctua- 
tions are small in both the clustered and gas phases — but 
increases significantly at T ~ 3 MeV where both phases 
coexists. 



IV. CONCLUSIONS AND OUTLOOK 

The quest for physical observables that are particu- 
larly sensitive to pasta formation remains elusive. In- 



O 



4 
3.5 

3 
2.5 

2 
1.5 

1 

0.5 






p/p =0.2 . 


■ Coexistence 

■ 1 


x p =o.3 ; 




Gas 

1 : 



2 4 6 8 10 12 14 16 18 20 
T [MeV] 



FIG. 14. The heat capacity of neutron-star matter as a func- 
tion of temperature and at a density of p/p = 0.2 and a 
proton fraction of a- p = 0.3. 



deed, even the existence of the pasta phase in the proton- 
deficient environment of the inner stellar crust remains 
an open question. In the present work we introduced 
a new model —the "Charged Ising Model" — to tackle 
some of these fundamental questions. The CIM is a 
lattice-gas model that while simple in its assumptions, 
retains the essence of Coulomb frustration that is re- 
quired to capture the subtle dynamics of the inner stel- 
lar crust. Monte Carlo simulations based on this model 
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are not as computationally demanding because the long- 
range Coulomb interaction — computed here exactly via 
an Ewald summation — was pre-computcd at the appro- 
priate lattice sites and then stored in memory prior to 
the start of the simulation. This represents an enormous 
advantage in simulating systems with a large number 
of particles as a function of density, temperature, and 
proton fraction. In this first — and mostly exploratory 
study — we were able to simulate systems with as many 
as S = 30 3 lattice sites as a function of temperature, den- 
sity, and proton fraction. Particular attention was placed 
on physical observables such as the pair correlation func- 
tion, the corresponding static structure factor, and the 
heat capacity; quantities that properly capture both den- 
sity and thermal fluctuations in the system. Note that 
the static structure factor displays a prominent coherent 
peak that occurs at a momentum transfer for which the 
probe (e.g., electrons in the case of protons and neutrinos 
in the case of neutrons) can most efficiently scatter from 
the density fluctuations in the system. The existence 
of pasta structures can therefore significantly reduce the 
electron or neutrino mean- free path in the stellar crust . A 
detailed study of electron and neutrino transport within 
the CIM will be forthcoming. We note that the very 
simple nuclear part of the Hamiltonian [Eqs. ([2]) and (J3J] 
essentially represents an Ising Hamiltonian for a spin-one 
system. That is, 

£ s 

Vkuclear = ~ X! ' ( 23 ) 

where e~ 16/3 MeV but now n,; = — 1, 0, 1 can take three 



values depending on whether the site is occupied by a 
neutron, is empty, or is occupied by a proton. To prop- 
erly simulate neutron-star matter, this Hamiltonian — 
together with the long-range Coulomb part — must be 
solved at constant magnetization rather than at constant 
magnetic field. Although much work has been done along 
the lines of the Ising model, the virtue of the CIM is that 
it incorporates the relatively simple spin-1 Ising model 
together with the challenging long-range Coulomb inter- 
actions. 

Finally, in the future we plan to use the CIM to carry 
out an analysis similar to the one recently reported in 
Ref. ■ In such a condensed- matter study a large tem- 
perature gap, i.e., kBT c < ksT<^ -Erjoui, was identified be- 
tween the melting temperature and the Coulomb energy 
where the system displays unconventional "pasta-like" 
behavior as a result of the strong frustration induced 
by the long-range interactions. Particularly relevant is 
the emergence of a pseudogap in the density of states 
that appears to mediate the transition from the Wigncr 
Crystal to the uniform Fermi liquid. We arc confident 
that the evolution of the pseudogap region as a function 
of proton fraction may help us prove the existence — or 
lack-thereof — of a pasta phase in the inner stellar crust. 
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